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\Q . Abstract 

For hyperbolic first-order systems of linear partial differential equations 
| (master equations), appearing in description of kinetic processes in physics, 

biology and chemistry we propose a new procedure to obtain their complete 
S closed-form non-stationary solutions. The methods used include the classical 

j> \ Laplace cascade method as well as its recent generalizations for systems with 

J^j ' more than 2 equations and more than 2 independent variables. As an example 

^ ■ we present the complete non-stationary solution (probability distribution) for 

Verhulst model driven by Markovian coloured dichotomous noise. 

Keywords: Master equations, hyperbolic systems, complete non-stationary 
solutions, kinetic processes, Verhulst model. 



1 Introduction 

This paper is devoted to a novel application of methods of explicit integration of 
hyperbolic linear systems of PDEs recently developed in [181 EE EQ] to an important 
class of dynamical nonlinear systems driven by a coloured noise. 

*This paper was written with partial financial support from the RFBR grant 06-01-00814 and 
the DFG Research Unit 565 'Polyhedral Surfaces" (TU-Berlin). 
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Modelling dynamical systems, in which stochastic behaviour is involved, as a rule 
leads to nonlinear stochastic differential equations for the dynamical variable or sets 
of dynamical variables. The comprehensive statistical treatment of these variables 
may be formulated with the aid of Liouville stochastic equation for the probability 
distribution (see [8]). 

Example 1. In this paper we consider as the simplest example the following 
one-dimensional dynamical system 



where x(t) is the dynamical variable, p(x), q(x) are given functions of x, a(t) is 
the random function with known statistical characteristics. The model (OQ) arises 
in different applications (see for example [H [10] and bibliography therein). An 
important application of this model consists in study of noise-induced transitions in 
physics, chemistry and biology. The functions p(x), q(x) are often taken polynomial. 
For example, if we set p(x) = p\X + p 2 x 2 , q(x) = q 2 x 2 p\ > 0, p 2 < 0, \p 2 \ > 
q 2 > 0, then the equation (jTJ describes the population dynamics when resources 
(nutrition) fluctuate (Verhulst model). In the following we will assume a(t) to be 
binary (dichotomic) noise a(t) = ±1 with switching frequency 2v > 0. As one can 
show (see [TOED]), the averages W(x,t) = (W(x,t)) and Wi(x,t) = (a(t)W(x,t)) 
for the probability density W(x, t) in the space of possible trajectories x(t) of the 
ODE satisfy a system of the form (JTj) (also called "master equations"): 



We suppose that the initial condition W(x, 0) = W (x) for the probability distribu- 
tion is nonrandom. This implies that the initial condition for Wi(x, t) at t = is zero: 
Wi(x,0) = (a(O)VF(x, 0)) = (a(0))W (x) = 0. The probability distribution W(x,t) 
should be nonnegative and normalized for all t: W(x, t) > 0, /f^ W(x, t) dx = 1. 

Example 2. Let us consider the following dynamical system driven by two sta- 
tistically independent Markovian dichotomous noises a(t) and f3(t): 



x = p(x) + a(t)q(x) 



(1) 



W t + (p(x)W) x + (q(x)W 1 ) x = 0, 

(WOt + 2vW y + (p(x)W 1 ) x + (q(x)W) x = 0. 



(2) 



x = p(x) + a(t)q(x) + (3{t)g(x). 



(3) 



The averaged probability density W(x, t) = (W(x,t))\ a ^ satisfies the following sys- 



tem of master equations: 



' W t + (p(x)W) x + (q(x)W 1 ) x + (g(x)P) x = 0, 
(Wi)t + 2vWi + ( P (x)W 1 ) x + a 2 (q(x)W) x + (g(x)Q) x = 0, 
(P) t + 2fiP + (p(x)P) x + (q(x)Q) x + 15 2 (g(x)W) x = 0, 
(Q) t + 2{n + u)Q + (p(x)Q) x + a 2 (q(x)P) x + (3 2 (g(x)W 1 ) x = 0, 
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where the auxiliary functions Wi(x,t), P(x,t) and Q(x,t) are some averages over 
realizations of the noises a(t), f3(i). They play the same auxiliary role as the func- 
tion Wi(x,t) in the system (T2]). We suppose that a(t) = ±a, (3{t) = ±j3 for any 
t. The characteristic switching frequencies of these random noises are 2v and 2fi 
respectively. We again suppose that the initial condition W(x, 0) = W (x) is non- 
random, so the initial conditions for Wi(x, t), P(x,t) and Q(x,t) at t = are zeros. 
The probability distribution W(x,t) should be nonnegative and normalized for all 
t: W(x, t) > 0, JZo W{x, t) dx = 1. 

Example 3. We can also consider nonlinear dynamical systems of higher order: 

i = f(x,y) + a(t)q(x,y), 
V = d(x,y) + P(t)s(x,y), 

where /, g, d, s are given functions. We use the same conventions for the noises 
a(t), (3{t) and the averaged probability distributions as in Example 2. The master 
equations for the main average P(x, y, t) and auxiliary averages Pi (a;, y, t), Q(x, y, t) 
and Qi(x,y,t) become 

' Pt + (fP) x + (gPi) x + (dP) y + (sQ) y = o, 

(P 1 ) t + 2uP 1 + (fP 1 ) x + a 2 (gP) x + (dP 1 ) y + (sQ 1 ) y = 0, 
Q t + + (fQ) x + (gQ 1 ) x + (dQ) y + (3 2 (sP) y = 0, 
t (QO* + 2(// + v)Q x + (fQ^ + a 2 (gQ) x + (dQ 1 ) y + (3 2 (sP 1 ) y = 0, 

The Cauchy initial value problem is formulated in the same way as above: P(x, y, 0) = 
P (x,y) is nonrandom; Px(x,y,t), Q(x,y,t) and Qi(x,y,t) at t = are zeros; 
P(x,y,t) should be nonnegative and normalized for all t: P(x,y,t) > 0, 
SZoJZoP{x,y,t)dxdy=l. 

Numerous publications (see [HI 1X0] for the bibliography) are devoted to solution of 
the system ([2]) asymptotically for t —>■ oo, that is to stationary solutions. A number 
of important phenomena of noise-induced phase transitions with applications in 
physics, chemistry and biology were discovered in this stationary case. 

There are only a few publications dealing with non-stationary solutions of (j2j). 
We refer to [16] and the recent paper [3], where some (incomplete) exact solutions 
of the system ([2]) for particular forms of the functions p(x) and q(x) were obtained. 

As we show in this paper, some interesting non-stationary kinetic equations 
(master equations) for probability distributions allow complete explicit closed-form 
solution of the general Cauchy initial value problem. These complete solutions are 
obtainable through a modification of the classical Laplace cascade method (see e.g. 
[USE]). This method is applicable to hyperbolic systems with two first-order linear 
PDEs in the plane (as (j2J) above) or a single second-order linear PDE in the plane. 
A preliminary closed-form complete solution for was obtained by this method in 
[19J. In that paper a much more general method of explicit integration, applicable 
to arbitrary hyperbolic higher-order linear systems (or a single higher-order linear 
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PDE) in the plane was developed. Later another generalization was proposed in 
|20j, it gives closed-form complete solutions for some special class of second-order 
linear hyperbolic equations with more than two independent variables. 

We give a brief account of the classical Laplace method as well as its new gener- 
alizations in Section [2j Section [3] is devoted to a detailed study of the system (j2J) for 
the simplest case of polynomial coefficients p(x) = p\X + p 2 x 2 , q(x) = q 2 x 2 , pi > 0, 
P2 < 0, \p2\ > q2 > (Verhulst model). We show that for an infinite sequence 
of values of the switching frequency v = p%, v = 2pi, v = 3pi, . . . , the complete 
explicit solution of the Cauchy problem is obtainable by our methods. 

In the final Section 0] we discuss future prospects and possible applications of 
our methods to more complicated systems of type (TjJ, (Q. 



2 Explicit integration of hyperbolic systems 
2.1 Laplace cascade method 

We give here only a special form of this method suitable for our purpose, see (3J [5J 
El [19] for more details. 

Suppose we are given a 2 x 2 first-order linear system of PDEs 

fvA fan a 12 \fv 1 \ i (b n b 12 



/x \ / \ /y 



i — i 1 1 ,+ r 11 r 1 (7) 

v 2 1 \a 2 i a 22 \v 2 \o 2 i o 22 ' 1 - ■ ' 




with aij = dij(x,y), bij = bij(x,y). We will suppose hereafter that (JTj) is strictly 
hyperbolic, i.e. the eigenvalues Xi(x, y), X 2 (x,y) of the matrix ( real and 

distinct. Let pi = {jp 11 {x,y),p 12 {x,y)), p 2 = (p 2 i(x,y),p 22 (x J y)) be the correspond- 
ing left eigenvectors: J^kPikUkj = \Pij- Form the following first-order differential 
operators Xj = D x — XiD y (the characteristic vector fields) and change the ini- 
tial unknown functions Vi to new characteristic functions iij = J^kPik^k- Then 
XiUi = J2k( x iPik)vk + EfcPifc ((vk)x - Xi(v k ) y ) = J2k, s Pik (aks - Xidks) (v s ) y + 
J2k,sPikbksVs+J2k(XiPik)vk=J2 s Vs (EkPikb ks + {XiPisfj = J2kUka ik (x,y), so we ob- 
tain the following characteristic form of the system (JTJ): 

f XiUx = a u (x,y)u 1 + a 12 (x,y)u 2 , 
\ X 2 u 2 = a 21 (x,y)u 1 + a 22 (x,y)u 2 , 

The characteristic system (jSJ), equivalent to ((Tj), is determined uniquely up to oper- 
ator rescaling Xi y)X and gauge transformations Ui — > gt(x, y)ui. It is easy 
to check that the gauge transformations to not change the Laplace invariants of the 
system h = X 2 (a n ) - Xx(a 22 ) - XiX 2 ln(ai 2 ) - X X {P) + Pa u + a 12 a 21 + (a 22 + 
X2(lna;i2) + P)Q and k = ct\ 2 a 2 \. Here P(x,y) and Q(x,y) are the coefficients of 
the commutator 

[X 1} X 2 \ = X X X 2 - X 2 X 1 = P(x, y)X x + Q(x, y)X 2 . (9) 
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These invariants h(x, y) and k(x, y) are just the classical Laplace invariants (cf. 
[HO El EH]) of the second-order scalar equation, obtained after elimination of u 2 from 
(jHD- Rescaling transformations of Xi change the Laplace invariants multiplicatively: 
h -> 7172/1, fc -> 7i72&- 

Starting from (jSJ) one can obtain two different (inequivalent w.r.t. gauge trans- 
formations) second-order scalar equations, eliminating either U\ or u 2 from ([8]). This 
observation gives rise to the Laplace cascade method of integration of strictly hyper- 
bolic systems in characteristic form (jSJ): 

(£1) If k vanishes then either au or a 2 i vanishes so the system becomes trian- 
gular: 

I Xiui = an iti, Qr I Xiiti = an iti + ai2 it2, /-^ 

\ X 2 u 2 = a 2 i «i + a 2 2 «2) 1 X 2 u 2 = «22 Wa- 

if we perform an appropriate change of coordinates (x, y) — > (x, y) (NOTE: for this 
we have to solve first-order nonlinear ODEs dy/dx = —\i(x,y), cf. Appendix in [7]) 
one can suppose X\ = D%, X 2 = Dy and obtain the complete solution of (TTUT) in 
quadratures: if for example a 12 = 0, then 

u 1 = Y(y)exp(—fa n dx), 

u 2 = exp (- / a 22 dy) (X(x) + j Y(y) exp (/(an dx - a 22 dy)) dy) 

where X(x) and Y(y) are two arbitrary functions of the characteristic variables x, 
y respectively. 

(C 2 ) If k 0, transform the system into a second-order scalar equation elimi- 
nating u 2 from (jHD: from the first equation 

u 2 = (Xiux - anMi)/ai2, (12) 

substitute this expressions into the second equation obtaining Lu\ = X 2 -^(XiUi — 

aniti) — a 2 \U\ — ^(XxUi — anUi) = 0. Now, using the commutator relation 0, 

we can represent Lu\ as Lu\ = {XiX 2 + f5\X\ + (3 2 X 2 + /3 3 )tti = (X% + /3 2 )(X 2 + 
fli)ui — hu\ = 0. From this form we see that this equation is equivalent to another 
2x2 system 

j X 2 ui = -/3i«i + u 2 , / 1Q x 

[ Ai« 2 = /l«i - fJ 2 U 2 . 

This new system (we will call it Xi-trans formed system) has the same characteristic 
form (jHD with different coefficients in the right-hand side. It also has new Laplace 
invariants hm, km, and it turns out that km equals to the invariant h of the original 
system. So if we have k^ = h = 0, we solve this new system in quadratures and 
using the same differential substitution ([121 we obtain the complete solution of the 
original equation Lu = 0. 

(£3) If again km 7^ 0, apply this Xx-transformation several times, obtaining a 
sequence of 2 x 2 characteristic systems TC(2), ^(3), • • • If on any step we get k( m ) = 
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0, we solve the corresponding system in quadratures and, using the differential 
substitutions f[T2"j) . obtain the complete solution of the original system. Alternatively 
one may perform X 2 -transformations, eliminating u± instead of u 2 on step (£2). In 
fact this ^-transformation is a reverse of the X\ -transformation up to a gauge 
transformation (see [lj). So we have (infinite in general) chain of systems 

. . . <— ri{-2) <— rt(-i) <— n — > /t(i) — > /t( 2 ) —>•... (14) 
and the corresponding chain of their Laplace invariants 

...,fc(_ 3 ), fc(_ 2 ), fc(-l), fc, k (1) = h, fc( 2 ), fc( 3 ),... (15) 

We do not need to keep the invariants fyj) in ffT5]) since ku\ = hu^iy If on any step 
we have kfjf\ = then the chains ( fl4l) and ( Tl5l) can not be continued: the differential 
substitution ( TT2l) is not defined; precisely on this step the corresponding system (jSJ) 
is triangular and we can find its complete solution as well as the complete solution 
for any of the systems of the chain fHM . 

As one may prove (see e.g. [3]) if the chain tHM is finite in both directions (i.e. 
we have fc(jv) — 0, k[-K) = for some N > 0, K > 0) one may even obtain a 
quadrature-free expression for the general solution of the original system: 

u x = c F + Cl F' +...+c N FW+d G + d 1 G' + ... + d K+1 & K+1 \ 

„ ~ (16) 

u 2 = e F + e x F> +. . .+e N F( N h f G + fiG' + . . . + f K+l G^ K+1 \ 

with definite Ci(x,y), di(x,y), ei(x,y), fi(x,y) and two arbitrary functions F(x), 
G(y) of the characteristic variables. Vice versa: existence of (a priori not complete) 
solution of the form ( JT6l) with arbitrary functions F, G of characteristic variables 
implies fc( a ) = 0, fc(_ r ) = for some s < N, r < K. So minimal differential 
complexity of the answer (TIB"]) (the number of terms in it) is equal to the number of 
steps necessary to obtain vanishing Laplace invariants in the chains (JHJ), ()15p and 
consequently triangular systems. Complete proofs of these statement may be found 
in [U t. 2], P2 E] for the case X\ = D x , X 2 = D y , for the general case cf. [El p. 30] 
and P]. 

We give a detailed example of application of this method in Section [31 
There were some attempts to generalize Laplace transformations for higher-order 
systems or the number of independent variables larger than 2, both in the classical 
time P [HI [15] and in the last decade [21 HE] • As one can show, all of them essentially 
try to triangulize the given system in some sense. A general definition of "generalized 
factorization" (triangulation) comprising all known practical methods was given 
in [18]. Unfortunately the theoretical considerations of [18] did not provide any 
algorithmic way of establishing generalized factorizability of a given higher-order 
operator or a given higher-order system. Below we present other approach for search 
of "generalized factorizations" resulting in explicit complete solution of some classes 
of hyperbolic systems. 
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2.2 Generalized Laplace cascade method for nxn hyperbolic 
systems in the plane 

Any nxn first-order linear system 

n n 

(vi) x = a ik(x,y)(v k )y + Yb ik (x,y)v k (17) 

k=l k=l 

with strictly hyperbolic matrix (a^) (i.e. with real and distinct eigenvalues of this 
matrix) is equivalent to a system in characteristic form 

XiUi = Y Oi ik (x, y)u k , (18) 

k 

as a straightforward calculation similar to that in the beginning of Section 12.11 im- 
mediately shows. 

Our generalization of the Laplace transformations consists in the following. 

(£i) For a given nxn characteristic system ffTB"]) choose one of its equations with 
a non- vanishing off-diagonal coefficient ^ 0, find u k = (XiUi — J2 s ^k a is u s)/&ik 
and substitute this expression into all other equations of the system. We obtain one 
second-order equation 

X k { — {XiUi a is u s )) - a kpU p - —(XiUi - a is u s )=0 (19) 

Qik S + k P ^k Qik B + k 

and n — 2 first-order equations 

X 3 u i ~ a Js u s —{XiUi - Y ®isU s ) = (20) 

s^k aik s^k 

for j i, k. 

(£2) The second step consists in rewriting the system (TT9]) . (|20|) in the following 
form with slightly modified unknown functions Tij = Uj+pj(x, y)u iy j 7^ i, k, = u iy 
new coefficients fl pq (x,y) but the same characteristic operators X p : 

Xi{—{X k Ui - J2(3i S u s )) -YlPkpUp - —{X k Ui ~Y PisUs) =0, (21) 
aik B+k p^k aik S+k 

XjUj - PjsUs ~ —{X k Ui - Pis^s) =0, j ^ i, k. (22) 

s^k aik s^k 

As one can prove this is always possible in a unique way. 

(£3) Introducing u k = -^{X k u,i — J2 s =i k j3i S u s ) rewrite ( [211) . ( 1221) as the trans- 
formed characteristic system 

X-' t U k ^2p fikptLpi 

X k Ui = £ s ^fc (3 is u s + a ik u k , (23) 



The reason of doing such generalized Laplace transformation consists in the fact 
that after it (or after a chain of such transformations) one may obtain a triangular 
system f[T8"j) . solve it in quadratures and doing the inverse steps (£3), (£2), (£1) 
obtain the complete solution of the original n x n hyperbolic system. 

Cf. [19] for the details, an example and the proof of correctness of the step 

2.3 Explicit solution of equations with more than two inde- 
pendent variables 

The method described in this Section is based on an idea given by Ulisse Dini in 
1902 (cf. [20J for references and details). We will limit here to a simple example 
showing the idea of the method. 

Let us take the following equation: 

Lu = (D x D y + xb x D z - D z )u = 0. (24) 

It has three independent derivatives D x , D y , D z , so the Laplace method is not 
applicable. On the other hand its principal symbol splits into product of two first- 
order factors: £1^2 + ^1^3 = £i(£2 + #£3)- This is no longer a typical case for 
hyperbolic operators in dimension 3; we will use this special feature introducing 
two characteristic operators X\ = D x , X 2 = D y + xD z . We have again a nontrivial 
commutator [Xi, X 2 ] — D z — X 3 . The three operators X^ span the complete tangent 
space in every point (x, y, z). Using them one can represent the original second-order 
operator in one of two partially factorized forms: 

L = X2X1 — X3 = X1X2 — 2X3. 

Let us use the first one and transform the equation into a system of two first-order 
equations: 

* = °~ {*&=*„. <*> 

Here comes the difference with the classical case dim = 2: we can not express u as 
we did in (fT2l) . But we have another obvious possibility instead: cross-differentiating 
the left hand sides of (|25j) and using the obvious identity [Xi,X 3 ] = [D X ,D Z ] = 
we get XiX 2 v = D x (D y + xD z )v = X 3 v = D z v or = D x (D y + xD z )v — D z v = 

(b x b y + xb x b z )v = (b y + xb z )b x v = x 2 x lV . 

Since we have now another second-order equation which is "naively" factorizable 
we easily find its complete solution: 

v — J <p(x, xy — z) dx + ip(y, z) 

where and ip are two arbitrary functions of two variables each; they give the 
general solutions of the equations X 2 4> = 0, Xiip = 0. 
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Now we can find u: 

U = I ( v ^ x + (D y + xD z )v dz^J + 6(y), 
where an extra free function 9 of one variable appears as a result of integration in 

(USD. 

So we have seen that such Dini transformations (I25I) in some cases may produce a 
complete solution in explicit form for a non-trivial three-dimensional equation (1241) . 
This explicit solution can be used to solve initial value problems for (1241) . 



3 Verhulst model 

Here we describe in detail the procedure of solution for the system ([2]). 

The characteristic operators and left eigenvectors of this 2 x 2-system are simple: 
Xi = I) t - XiD x , Ai, 2 = -p(x) ± q(x), pu = p 2 i = P22 = 1, P12 = -1- The 
characteristic system ([8]) for the new characteristic functions U\ = W — W\, u 2 = 
W + W x is 

| X]U\ = -{p x - q x + v)u\ + uu 2 , ^ 
\ X 2 u 2 = uui- (p x + q x + v)u 2 . 

The Laplace invariants are h = u 2 — [p X xq 2 (p + <?) +plq 2 —Pxq x q{?>p + q) — q xx pq{p + 
q) — q x p(2p + q)]/ q 2 , k = z/ 2 , so if u, p(x) and q(x) satisfy a second-order differential 
relation h = 0, one can solve ([2]) in quadratures. Especially simple formulas may 
be obtained for polynomial p(x) = p\x +p 2 x 2 , q(x) = q 2 x 2 : in this case k = u 2 , 
h = h(-2) = v 2 — p\ so if v = pi, one may solve ([2]) explicitly. It is convenient at 
this point to use the dimensionless variable r = ut; so we have to change t i-> t, 
v i — ► 1, p\ i — ► 1 and change p 2 , q 2 respectively. For simplicity we will still use the 
same notations p 2 , q 2 . 

After the necessary transformation, described in Section 12.21 we obtain the fol- 
lowing quadrature-free expression for the complete solution of the system (J2J): 

W = % [F'{x) - F(x) + G'(y) - G(y)] , 
x z 

Wi = ^ {-q 2 xG\y) + (1 + p 2 x)G{y) + q 2 xF\x) ( 2T ) 
+ (l+p 2 x)F(x)}, 

where x = — t + In , , , -. — t-, y = — t + In . . , x — r- are the characteristic variables 

(X 2 x = 0, X\y = 0) and F, G are two arbitrary functions of the corresponding 
characteristic variables. 

For the case v 2 ^ p\ we can compute other Laplace invariants of the chain f|T5l) : 
h(i) = = v 2 - 4pf, /i( 2) = /i(_4) = v 2 - 9pj, h( 3 ) = /i(_ 5) = v 2 - 16pf , etc., so 
for the fixed p[x) = p\X + p 2 x 2 , q(x) = q 2 x 2 and v = ±pi, v = ±2p 1; v = ±3pi, 
. . . one can obtain closed- form quadrature-free complete solution of the system ([2]), 
with increasing complexity of the answer ( TL6]) . 
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Now we demonstrate how the formulas (127)) may be used to solve the Cauchy 
initial value problem. For this set r = inside the variables x, y and equate 



l+(P2+<?2) 



Hi V = ln T- 



{v2-q2)x 



one 



W(x,0) = Wo(x), Wi(x, 0) = 0. Since now x = In 
can express the derivatives F' = dF/dx, G' = dG/dy as F' = G' — and 

obtain from (127)) a system of two linear ODEs for F(x) = F(x(x)), G(x) = G(y(x)). 
It may be solved explicitly (see an explanation of this fact in Section HJ) for any 
W (x): 



F(x) 
G{x) 



i 



2ql(l + { P2 +q2)x) 

1 

2<?l(l+(P2-g 2 )z) 



x W (6>) 

co e 



d0 + (l + q 2 x)J*W o {6) d9 



xj: ^de + ( q2 x-i)j: i w (e)de 



(28) 



Perform now the inverse substitution F(x) = F(x(x)), G(y) = G(x(y)) (for r = 0) 
to find the "true" functions F(x), G(y) suitable for substitution into (|2~7j) for any r. 
This final form of the explicit solution of the Cauchy problem is: 



W(x,t) 



1 



2q 2 x : 



[h(y) - h(x)} + 



W {y) 



+ 



2(e r (l + (p 2 - q 2 )x) - x(p 2 - q 2 )) 2 
W (x) 



2(e r (l + (p 2 + q 2 )x) - x(p 2 + q 2 )) 2 



W 1 (x,t) 



h{y)-h{x) 

2q 2 x 3 
W (y) 



[e T - l)p 2 x - 1 



+ 



2q 2 x 2 
W (x 



[h{v)-h{x)\ 



2(e T (l + (p 2 - q 2 )x) - x(p 2 - q 2 )) 2 2(e r (l + (p 2 + q 2 )x) - x(p 2 + q 2 )) 2 ' 

where x = {eT{1+{p2+q2 x )x) _ x{p2+q2)) , y = ^{i+(p 2 - q2 X )x)-x{p 2 -q2))' Ix ^ z ' = & w o(°) 
h(z) = fc ^j^-dO, Co and c\ may be chosen arbitrary 
One can check that x < y for all t > 0, x > 0. 

We get an especially simple form of this solution for the initial distribution 
Wq(x) = 5(x — x*) with some fixed initial state x(0) = x* > 0: 



W{x,t) 



5(x — x*) 



+ 



6(y - x*) 



2(e T (l + (p 2 + q 2 )x) - x(p 2 + q 2 )) 2 2(e T (l + (p 2 - q 2 )x) - x(p 2 - q 2 )) 



+ 



H(y - x*) - H(x - x* 
2q 2 x 2 



Here H(z) = f*^ 5(9) d6 is the Heaviside function. 

According to the standard formula <5(0(x)) = 5(</> _1 (0))/<//( 

Sjx - X*) _ &{ x ~ l-(p 2 +g 2 KV-l)^;) 

2(e T (l + (p 2 + q 2 )x) - x{p 2 + q 2 )) 2 " 2e T 



(0)) one gets 
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5{V- X *) 5 ( X ~ l-(p2-q2)fr-l)xj 

— - — , so we see that this 



2(e T (l + (p 2 -q 2 )x) -x{p 2 -q 2 )) 2 2e T 
simple solution (and consequently the complete solution) obviously obeys the neces- 
sary physical requirements of positivity and normalization: W(x, t) > 0, W(x, t) dx 
1. Asymptotically, for r — > oo, this solution exponentially fast converges to the sta- 
tionary probability distribution W OQ (x) = outside the interval | P2 l g2 < x < \ P2 + Q2 \ 
and Woo (a;) = l/{2q 2 x 2 ) inside this interval. 



4 Concluding remarks and future prospects 

There is an algorithmic possibility to obtain closed-form solutions of the Cauchy 
problem for the more complicated cases v = mpi, m = 2,3,4, .. . in the Verhulst 
model. The respective classical form ffl6|) is algorithmically obtainable with the 
methods of Section I2TT1 Since the orders of derivations of F(x), G(y) in the right- 
hand sides of f[T6"j) are proportional to the integer coefficient m in the relation v = 
mpi, directly assigning W(x, 0) = W (x), Wi(x, 0) = in this formula for r = will 
result in a linear system of ODEs for Fix), G(y) of high order with nonconstant 
coefficients. Much more efficient is to use the transformations (fl2l) directly: simply 
recalculate the Cauchy data for the new functions u 2 on step (£2), using ( T131) . until 
we get (after m steps) a triangular system, solve this system for the recalculated 
Cauchy data and then use the inverse ^-transformations to get the solution of the 
original system. This also explains why we could find the solution (128]) in the case 
v = pi in Section O 

Methods, described in Sections 12. 2} 12.3} suggest that one can also investigate 
systems 01]), and classify completely integrable cases for special forms of their 
coefficients p(x), g(x), f(x,y), q(x,y), d(x,y), s(x,y) and switching frequencies /x, 
^.Systematic investigation of such integrable cases will be reported in subsequent 
publications. 
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